
rm(list=ls())

REP.index = 1 

source("2.0.0 MyFunc.R")

load("100716MSCM.RData") # Change accordingly
data = data[,2:ncol(data)]  # delete intercept
pL0  = pL0 - 1
source("3.2 PMLE.R")
source("3.5 CallDR.R")
REP.each = 5

max.steps = c(50)
btstrap.name  = c(paste("btstrap=",1:REP.each,sep=""),"point.est")
para.name = c("theta", "phi", "gop", paste("eta",c("0*","1*")),
              "nll time", "contrast")
results  = array(dim=c(10,REP.each + 1,length(max.steps), 7,pL0+3), dimnames = 
                    list(1:10,btstrap.name, max.steps, para.name,
                         c("*=0", "*=1", "*=.", 1:pL0)))
names(dimnames(results)) = c("Simulation number","Bootstrap number",
                             "Max steps","Param","subscript")

results.dr.MLE = results.dr.logit = results.pmle = results


for(i in 1:length(max.steps)){
    # i= 1; seed.start = 1# DEBUG
    append.to.result = function(results, est, btstrap.arg, seed.start){
        results[seed.start,btstrap.arg,i,1, ]    = est$theta.coef
        results[seed.start,btstrap.arg,i,2, ]    = est$phi.coef
        results[seed.start,btstrap.arg,i,3, ]    = est$gop.coef
        results[seed.start,btstrap.arg,i,4:5, ]  = est$eta.coef
        results[seed.start,btstrap.arg,i,6, 1:3] = c(est$nll, est$time, est$n01)
        results[seed.start,btstrap.arg,i,7, 1:3] = est$contrast
        return(results)
    }
    min.nll = 1000; idx.min.nll = 1
    point.mle.lst = list()
    for(seed.start in 1:10){
        point.mle = PMLE(data, count, K, seed = seed.start, max.step = max.steps[i], 
                         thres = 1e-4) 
        point.mle.lst[[seed.start]] = point.mle
        results.pmle = append.to.result(results.pmle, point.mle, btstrap.arg = "point.est",seed.start=seed.start)
        if (point.mle$nll < min.nll){
            idx.min.nll = seed.start
            min.nll = point.mle$nll
            point.MLE = point.mle
        }
    }
    print(paste0("Minimum nll is: ", min.nll, ";It is the ",idx.min.nll,"th iteration"))
    
    for(seed.start in 1:10){
        point.dr.est.MLE = DREst(data, cnt= count, K=K, seed=seed.start, 
                                  max.step = max.steps[i], thres = 1e-4, 
                                  mle=point.MLE, outcome.type = "MLE", dr.warm = NULL)
        results.dr.MLE = append.to.result(results.dr.MLE, point.dr.est.MLE, btstrap.arg  = "point.est",seed.start=seed.start)
        if (point.dr.est.MLE$nll < 0.005){ # Here point.dr.est.MLE$nll is the DR objective value
            # Copy and paste this result and stop
            if (seed.start < 10){
                for (j in (seed.start+1):10){
                    results.dr.MLE = append.to.result(results.dr.MLE, point.dr.est.MLE, btstrap.arg  = "point.est",seed.start=j)
                }
            }
            break
        }
    }
    print(paste0("y.hat.type = MLE, Minimum DR Value is: ", point.dr.est.MLE$nll,"; It is the ",seed.start,"th iteration"))
    
    for(seed.start in 1:10){
        point.dr.est.logit = DREst(data, cnt= count, K=K, seed=seed.start, 
                                    max.step = max.steps[i], thres = 1e-4, 
                                    mle=point.MLE, outcome.type = "logit", dr.warm = NULL)
        results.dr.logit = append.to.result(results.dr.logit, point.dr.est.logit, btstrap.arg = "point.est",seed.start=seed.start)
        if (point.dr.est.logit$nll < 0.005){ # Here point.dr.est.MLE$nll is the DR objective value
            # Copy and paste this result and stop
            if (seed.start < 10){
                for (j in (seed.start+1):10){
                    results.dr.logit = append.to.result(results.dr.logit, point.dr.est.logit, btstrap.arg = "point.est",seed.start=j)
                }
            }
            break
        }
    }
    print(paste0("y.hat.type = logit, Minimum DR Value is: ", point.dr.est.logit$nll, "; It is the ",seed.start,"th iteration"))
    
    ######
    # Start bootstrapping 
    for (seed in 1:REP.each){
        # seed = 1 # DEBUG
        min.nll = 1000; idx.min.nll = 1
        for(seed.start in 1:10){
            SEED  = seed + (REP.index-1)*REP.each
            set.seed(SEED)
            space = rep(1:length(count),count)
            index = sample(space,size=length(space),replace=TRUE)
            count.btstrap = as.numeric(table(factor(index,levels=
                                                        1:length(count))))
            
            if (seed.start == 1){
                mle.warm.i = point.MLE
            }else{
                mle.warm.i = point.mle.lst[[seed.start]]
            }
            mle = PMLE(data, count.btstrap, K, SEED + seed.start + (seed-1) * REP.each,
                       max.step = max.steps[i], thres = 1e-4, mle.warm=mle.warm.i)

            results.pmle = append.to.result(results.pmle, mle, btstrap.arg = seed, seed.start = seed.start)
            
            if (mle$nll < min.nll){
                idx.min.nll = seed.start
                min.nll = mle$nll
                MLE = mle
            }
        }
        print(paste0("btstrap number: ", seed, "; Minimum nll is: ", min.nll, "; It is the ",idx.min.nll,"th iteration"))
        
        for(seed.start in 1:10){
            SEED  = seed + (REP.index-1)*REP.each
            set.seed(SEED)
            space = rep(1:length(count),count)
            index = sample(space,size=length(space),replace=TRUE)
            count.btstrap = as.numeric(table(factor(index,levels=
                                                        1:length(count))))
            if (seed.start == 1){
                dr.warm = point.dr.est.MLE
            }else if (seed.start == 2){
                dr.warm = MLE
            }else{
                dr.warm = NULL
            }
            dr.est.MLE = DREst(data, cnt= count.btstrap, K=K, seed=SEED + seed.start + (seed-1) * REP.each, 
                                      max.step = max.steps[i], thres = 1e-4, 
                                      mle=MLE, outcome.type = "MLE", dr.warm = dr.warm)
            results.dr.MLE = append.to.result(results.dr.MLE, dr.est.MLE, btstrap.arg  = seed,seed.start=seed.start)
            if (dr.est.MLE$nll < 0.005){ # Here point.dr.est.MLE$nll is the DR objective value
                # Copy and paste this result and stop
                if (seed.start < 10){
                    for (j in (seed.start+1):10){
                        results.dr.MLE = append.to.result(results.dr.MLE, dr.est.MLE, btstrap.arg=seed, seed.start=j)
                    }
                }
                break
            }
        }
        print(paste0("btstrap number: ", seed, ";y.hat.type = MLE, Minimum DR Value is: ",dr.est.MLE$nll,"; It is the ",seed.start,"th iteration"))
        
        for(seed.start in 1:10){
            SEED  = seed + (REP.index-1)*REP.each
            set.seed(SEED)
            space = rep(1:length(count),count)
            index = sample(space,size=length(space),replace=TRUE)
            count.btstrap = as.numeric(table(factor(index,levels=
                                                        1:length(count))))
            if (seed.start == 1){
                dr.warm = point.dr.est.logit
            }else if (seed.start == 2){
                dr.warm = MLE
            }else{
                dr.warm=NULL
            }
            dr.est.logit = DREst(data, cnt= count.btstrap, K=K, seed=SEED + seed.start + (seed-1) * REP.each, 
                                max.step = max.steps[i], thres = 1e-4, 
                                mle=MLE, outcome.type = "logit", dr.warm = dr.warm)
            results.dr.logit = append.to.result(results.dr.logit, dr.est.logit, btstrap.arg  = seed,seed.start=seed.start)
            if (dr.est.logit$nll < 0.005){ # Here point.dr.est.logit$nll is the DR objective value
                # Copy and paste this result and stop
                if (seed.start < 10){
                    for (j in (seed.start+1):10){
                        results.dr.logit = append.to.result(results.dr.logit, dr.est.logit, btstrap.arg=seed, seed.start=j)
                    }
                }
                break
            }
        }
        print(paste0("btstrap number: ", seed, ";y.hat.type = logit, Minimum DR Value is: ",dr.est.logit$nll,"; It is the ",seed.start,"th iteration"))
    }
}   

save(results, results.pmle,results.dr.MLE,results.dr.logit,
     file=paste("./RDataFiles/072221_btstrap_results_SEED",REP.index,".RData",sep=""))










